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Abstract 

Loop invariants play a very important role in proving correctness of programs. In this paper, we 
address the problem of generating invariants of polynomial loop programs. We present a new 
approach, for generating polynomial equation invariants of polynomial loop programs through 
computing vanishing ideals of sample points. We apply rational function interpolation, based 
on early termination technique, to generate invariants of loop programs with symbolic initial 
values. Our approach avoids first-order quantifier elimination and cylindrical algebraic decompo- 
sition(CAD). An algorithm for generating polynomial invariants is proposed and some examples 
are given to illustrate the algorithm. Furthermore, we demonstrate on a set of loop programs 
with symbolic initial values that our algorithm can yield polynomial invariants with degrees 
high up to 15. 

Key words: Program Verification, Polynomial Loop Programs, Invariant Generation, 
Vanishing Ideals 



1. Introduction 

Loop invariant generation plays a central role in program verification. An invariant 
of a loop program at a location is an assertion over the program variables that is true of 
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any program state reaching the location. Loop invariants are helpful for program analysis 
and verification. 

Since the late seventie s of the 20th century, many methods have been proposed to gen- 
erate loop i nvariants. In (IGerman and Wegbreitj . ll975tlKatz and Manna . 119761: Wcgbrcit. 



19741 119751 ). difference equation solving techniques were used to generate loop invariants. 



However, this technique i s diffi c ult to apply in general, since diff e rence equations are gen- 
erally hard to solve. In (jKarrl . 119761: ICousot and Cousotl Il977t ICousot and Halbwachsi 
19781 ). abstract interpretation techniques were applied to finding linear equation or in- 



equality invariants. 

Based on some previous work, Mullcr-Olm and Seidl ( Miiller-Olm and Seidl . 2004al lbl) 
generated polynomial equation invariants of loop programs with affinc assignments using 
linear algebra techniques. 

Recently, the constraint-based methods become dominant in invariant generation. 
These methods require to preset a template of the invariant as a polynomial equation or 
inequality with unknown coefficients, and the initiation and consecution conditions for 
the invariants generate constraints on the unknown coeffi cients. Then a solution to the 
constraint system yields invariants. In ( Colon et al. . 20031) . Farkas' Lemma was applied 



to ge nerating linear inequality invariants using non-linear constraint solving. In ( Kapurl . 
20041 ). Kapur proposed an ap proach based on quantifier elim ination to generate polyno- 



mial equation invariants. In ( Sankaranaravanan et al. . 120041) . the polynomial-scale con- 



secution of inductive invariants was first defined, and the polynomial equation invariants 
satisfying poly nomial-scale consec ution were computed using an extended Grobner ba- 
sis algorithm. ( Rebiha et all [2008h proposed a complete method using multi-parametric 
constraints to generate polynomial invariants that satisfy polynomial-scale consecution. 
To generate polynomial equation or inequality invaria nts of loop programs with guard 
conditions and branches, Chen et al. ( Chen et al. . 20071 ) applied the techniques of solving 
semi-algebraic systems. 

Rodrigucz-Carbonell and Kapur ( Rodriguez-Carbonell and Kapur . 2004 . 20071) first 
proved that the polynomial equation invariants have the algebraic structure of an ideal, 
and they proposed a fixpoint procedure for find ing all polynomial equation invariants us- 
ing Grobner bases a nd quantifier elimination. ( Kovacsl 2007 : Kauers and Zimmermann , 



2008t iKovacsl 120081) proposed complete algorithms to generate polynomial equation in- 



variants for a restricted class of linear (P-solvable) loops. 

In this paper, by computing vanishing ideals of program sample points, we present a 
new method for generating polynomial invariants of polynomial loop programs in which 
the guard conditions and assignments are polynomials in the program variables. 

Recall that a multivariate polynomial in n variables with total degree bound e has 
at most (™ +e ) distinct terms. Therefore, to compute the invariants with a given degree 
bound e, we first get no more than (™^ e ) sample points by executing the loop program, 
where n is the number of program variables. Then we apply Buchberger-Mollcr algorithm 
to compute the vanishing ideal of these sample points as candidate invariants (a candi- 
date may not be real invariant). Subsequently, the problem of verifying the candidate 
invariants can be translated into that of determining divisibility between multivariate 
polynomials, and a practical probabilistic method is presented to exclude non-invariants 
quickly. Finally, we can either generate the polynomial invariants or conclude that the 
polynomial invariants with degree < e' do not exist, where e'(< e) is the minimal de- 
gree of the polynomials in the vanishing ideal. Moreover, rational function interpolation 
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method, combining variab le by variable interpolation with early termination technique 
(IKaltofen and YansJ . 120011) is applied to generating invariants of loop programs with 
symbolic initial values. 

The rest of the paper is organized as follows. In Section 2, we recall the notions of 
vanishing ideals for finitely many points, transition systems and (inductive) invariants. 
In Section 3, we present an efficient method to generate polynomial equation invariants 
for polynomial loop programs with initial values, and in Section 4, an algorithm and 
some examples are given. We conclude our results in Section 5. 



2. Notation and Definitions 

2.1. Vanishing Ideals of Finitely Many Points 

This section contains a collection of definitions and facts about vanishing ideals of 
finitely many points. 

Throughout this paper, let K. be a (commutative) field of characteristic zero , K\xi, . . . , x n ] 
be the ring of polynomials in n indeterminates Xi,...,x n over K., the term order a in 
K,[xi, . . . , x n ] be the graded lexicographic order, and deg(/) denote the total degree of a 
polynomial / £ K\x\, . . . , x n }. 

Definition 1 (Ideal of Polynomials). A set X C K\xi, . . . , x n ] is an ideal in K\xi, . . . , x n ] 
if for any f,g£X and h 6 K\x\ , . . . , x n ] , we have / + g E X and / ■ h 6l, 

For hi, . . . ,h r € JC[xi, . . . , x n ], we denote by {hi, . . . , h r ) the smallest ideal contain- 
ing hi,...,h r , i.e. 



(hi,...,h r ) = < y~]fjhj | fi £ K[xi,.. 



k t=i ) 

If X = (hi, . . . , h r ), we say that X is an ideal generated by hi, . . . , h r and that hi, . . . , h r 
is a basis of X. 

By Hubert' Basis The orem, any ideal in K\xi , . . . , x n ] has a finite basis. The classical 
Buchberger's algorithm ( Buchberger , 19851) can be applied to computing Grobner bases 
of ideals. 

Definition 2 (Vanishing Ideal of Finitely Many Points). Let A be a finite subset of K, n . 
The vanishing ideal of the point set A is the ideal 

X(A) = {f &K[xi,...,x n ] | f{a)=0, for all a G A} 

of all the polynomials that vanish on each point in A. 

Buchberger and Moller presented an algorithm (jMoller and Buchberger , Il982h to com- 
pute the reduced Grobner bases of the vanishing ideals for finitely many points, via 
Gaussian elimination on a generalized Vandermonde matrix. 



Remark 1. As stated in (jMarinari et all 119931) . Buchberger- Moller algorithm is of poly- 
nomial time complexity 0(n 2 S*), where n is the dimension of the affine space JC n and s 
is the number of points. 
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2.2. Transition Systems and Invariants 

The classical way to represent programs is the use of transition systems. 

Definition 3 (Transition System). A transition system T is a tuple {V, L, T, Zo, O), where 

• V — {x\, ... , x n } is a finite set of program state variables; 

• L is a set of locations; 

• 1~ is a set of transitions, where each transition r £ T is a tuple of the form 

\h,h,Pri 9t), 

such that 

• h,l2 & L are the pre- and post- locations of r, respectively; 

• p T is a transition relation, i.e., a first-order formula, over V U V, where the prime 
version V' := {x[ , . . . , x' n } of V represents the next-state variables. Here x\ is the new 
variable introduced to stand for the value of Xi after the assignment. For example, 
the assignment x := x + 1 can be written as x' — x + 1; 

• g T is the guard condition of the transition r or of p T . Only if g T holds, the transition 
can take place. 

• Location Iq G L is an initial location, and the initial condition is a first-order formula 
over V. 

As an example, a loop program 

where O 
Z : while guard do 

xi := Pi(xi,...,x n ); 

Xn . Pn(xi : ■ • ■ , X n *) : 

end while 

can be translated easily into a transition system 

( ..,*„}, {Z}, {r}, {Z}, 6 ) 

where 

r = (l,l,x[ = Pi(xi,. . .,x n ) A •■• Ax' n = P„(xi,. . . ,x n ), guard). (2) 

A polynomial loop program is a loop program (1) where all the Pi are polynomials in 
X\, . . . ,x n , and the guard is represented by a conjunction of polynomial inequalities. 

By an assertion, we mean a first-order formula over the program variables. A state 
of a transition system is an interpretation of the program variables as values from the 
corresponding domains. We use the notation s (= ip to denote that a state s satisfies an 
assertion ip. We will also write tpi \= (f2 for two assertions ip±, tp2 to represent that if2 is 
true at least in all the states in which ipi is true. 

Next, we introduce the notions of (inductive) invariants for transition systems. 
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Definition 4 (Invariant). Let T = (V, L, T, lo, @) be a transition system. An invariant 
of the system T at location I <E L is an assertion over V which holds at all reachable 
states at location I. An invariant of the system T is an assertion over V that holds at 
all locations. 

Definition 5 (Inductive Invariant). Let T = (V, L, T, @) be a transition system and D 
the domain of assertions. An assertion map for T is a map rj : L — > D that associates 
each location of the transition system with an assertion. We say that rj is inductive if the 
Initiation and Consecution conditions hold: 
Initiation: |= 77(^0); 

Consecution: For each transition r = (Zi, I2, p T > 9t)i we have 

n(h)(V) A pr{V, V) A g T (V) |= r,{l 2 ){V), (3) 

where T)(l2)(V) represents the assertion 77(^2) with the current state variables x\, . . . ,x n 
replaced by the next state variables x'i, . . . , x' n , respectively. 



It is a well-known result ([Flovdl 119671 ) that if rj is an inductive assertion map then rj{l) 



is an invariant at location I for each / G L. 

Remark 2. When no confusion arises, we write the condition (3) simply as 

V (V)Ap T (V,V')Ag T (V)h V (V). 

In this paper, we are interested in finding inductive invariants of the form p(V) = 0, 
where p(V) is a polynomial in the program variables. For brevity, we shall use r)(V) to 
denote both the assertion p(V) = and the polynomial p(V). 

In the sequel, we will use the following strong er but more practi cal consecution con- 
dition defined in ( Sankaranaravanan et al. . 2004 , Definition 13) or ( Rebiha et al. , 20081 
Definition 4). 

Definition 6 (Polynomial-Scale Consecution). Let r = (hjfo, Pt:9t) be a transition 
and r\ be an assertion map. We say that r\ satisfies polynomial- scale consecution for t if 
there exists a polynomial q(V) such that 

Pr H (v(V) - q(V) ■ ri(V) = 0). 

In particular, if deg(<7) = 0, polynomial-scale consecution reduces to constant-scale con- 
secution. 

From Definition 6, to verify whether rj satisfies polynomial-scale consecution, it suffices 
to check whether r)(V) can be divided by ^(V), i.e., n{V) | r)(V). 



3. Generating Invariants of Polynomial Loops with Initial Values 

In this section, we present an approach to generate polynomial equation invariants for 
polynomial loop programs with initial values. 
Our idea is as follows. 

Step 1: We run the loop program and get a set S of sample points by recording the 
values of system variables at each location. 



5 



Step 2: We apply Buchberger-Moller algorithm to compute a Grobner basis {771, . . . , rj r } 
of the vanishing ideal 1(S) of S, and take the polynomials in the basis of I(S) (or more 
exactly, the corresponding polynomial equations) as candidate invariants. 

Step 3: For each candidate invariant, we determine whether it is an invariant of the 
given program by checking polynomial-scale consecution. 

3.1. Determining Divisibility between Multivariate Polynomials 

We now describe Step 3, the key part of our idea in details, that is, how to verify effi- 
ciently whether a polynomial satisfies polynomial-scale consecution. Clearly, this problem 
is equivalent to checking divisibility between multivariate polynomials. 

To fulfill this task, a straightforward way is to apply directly the multivariate polyno- 
mial division algorithm. In order to find out all polynomial invariants among the ideal 
basis 771,..., r) r , one will need to apply multivariate polynomial division algorithm r 
times. Actually in our experiments, we find that only a few of polynomials in the ideal 
basis are invariants, i.e., the number of invariants is much less than r. For example, in 
Example 1 to be presented, there are six candidate invariants while only one of them 
is actually an invariant of the given program. Taking this special case into account, in- 
stead of directly applying multivariate polynomial division algorithm, we present a high 
probability algorithm to determine multivariate polynomial divisibility, combining linear 
transformation and univariate polynomial division. The merit of this technique lies in 
reducing the computational complexity. 

The idea of our high probability method for determining multivariate polynomial 
divisibility is based on the following theorem. 

Theorem 1. Let f,g £ K,[xi, . . . ,x n ] with deg(/) = e\ < dcg(g) = e 2 , and e\ > 0. 
Suppose that B 2 , ■ • ■ , B n ,p2, . . . ,p n are chosen randomly and uniformly from a finite 
subset W C K. Let f,g£ IC[Z] be the univariate polynomials (in a new indeterminate Z) 
constructed from / and g, respectively, as follows 

f(Z) =f(Z,B 2 Z-p 2 ,...,B n Z-p n ) 
g(Z) =g(Z, B 2 Z -p 2 ,..., B n Z - p n ). 

If f\g, then f\g; or equivalcntly, if / \ g, then / \ g. 

Proof. If f\g, then there exists a nonzero polynomial h £ /C[a;i, . . . , x n ] such that g = f-h. 
From (4), we get g(Z) = f(Z) ■ ~h(Z) where 

h(Z) = h(Z, B 2 Z -p 2 ,..., B n Z - p n ). 

So we have f\g. □ 

Suppose that {771, . . . ,rj r } is a set of candidate invariants. By Definition 6, to check 
whether r\i satisfies polynomial-scale consecution it suffices to determine if r/i(V)\rii(V') 
for each i. Denote by f]i{Z) and fji(Z') the univariate polynomials as constructed in (4) 
from r)i(V) and r]i(V), respectively. According to Theorem 1, we have the following 
observations: 

(1) If fji(Z) \ fji(Z') then rji(V) \ rji(V'), i.e., r]i does not satisfy polynomial-scale conse- 
cution and then is not an invariant of the program. In such a way, the polynomials 
that are not invariants can be removed quickly from the invariant candidates by 
checking univariate polynomial divisibility. 




(» 



(2) It is not necessarily true that if fjj(Z)\fjj(Z') then T]j(V)\r]j(V'). Therefore, in 
the case fjj(Z)\fjj(Z'), to check whether r)j is really an invariant we need to ap- 
ply multivariate polynomial division to verify if rjj(V)\rjj(V'). However for fjj(Z) 
and fjj(Z') as constructed in (4), we will show in Theorem 2 that if rjj(V) f?7j(V') 
then rjj (Z) \ i)j (Z') with high probability. In other words, if the multivariate polyno- 
mial rjj does not satisfy polynomial-scale consecution, then with very low probabil- 
ity the corresponding univariate polynomials fjj{Z) and fjj(Z') satisfy fjj(Z)\fjj(Z'). 
Before presenting the main theorem, we need several lemmas. 



Lemma 1 (Schwartz-Zippel Lemma. ISchwarta (|1980f )). Let / £ K[xi, . . . ,x n ] be a non- 



zero polynomial with deg(/) = e > 0. Let If be a finite subset of JC of cardinality \W\ 
and let r-y, r%, . . . , r n be chosen randomly from W. Then we have the following probability 
estimate: 

Prob(/(ri,ra,...,r n )=0) < ^- 
Proof. In the univariate case, the proof follows easily from the fact that a univariate 



polyn omial of degree e has no more than e roots. The reader can refer to ISchwartz 
(jl980h :br the proof in the multivariate case. □ 



The probability estimate in Lemma 1 will be needed to show that the multivariate 
polynomial and its associated univariate polynomial as constructed in (4) have the same 
degree with high probability. 

Lemma 2. Let / € JC[xx, . . . ,x n ] with deg(/) = e > 0, W be a, finite subset of JC 
of cardinality \W\, and P>2, . ■ . , B n be distinct points chosen randomly from W. Let 
/ = f(Z, B2Z, . . . , B n Z) 6 K[Z]. Then we have the following probability estimate: 

Prob(deg/ = e)>l--|-. 

Proof. The polynomial / can be partitioned as two parts / = /1 + $2, where f\ consists 
of the terms in / with degree = deg(/) and fi consists of the terms in / with degree < 
deg(/). If deg(/) < deg(/), then we have h(Z, B 2 Z, B n Z) = 0. Then the probability 
estimate follows from Lemma 1. □ 

The following lemma gives an estimate of the probability that two univariate polyno- 
mials constructed as in (4) from two coprime multivariate polynomials remain coprimc. 

Lemma 3. Let f,g 6 K[x\, . . . , x n ] with deg(/) = e\, deg(g) = e2 and gcd(/, 5) = 1. 
Let W be a, finite subset of JC of cardinality \W\. Suppose that B2, ■ ■ ■ ,B n £ JC are given 
such that 

deg(f(Z,B 2 Z,...,B n Z)) = ei and dcg(g{Z, B 2 Z, . . . , B n Z)) = e 2 

where Z is a new indeterminate. Let P2, ■ ■ ■ ,p n be n — 1 distinct points chosen randomly 
and uniformly from W. As in (4), let 

f{Z) = f{Z, B 2 Z -p 2 ,..., B n Z - p n ) and g(Z) = g(Z, B 2 Z -p 2 ,..., B n Z - p n ). 

Then 

Prob(gcd(/(Z),fl(2)) = l)>l-^. 
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Proof. The proof is similar to that of Lemma 2.2 in (jKaltofen and Yand . l2007l) . For new 
variables Z, 012, ■ ■ ■ , a n we define the map: 

4>: JC[xi,x 2 , ...,x n ] -t K.[Z,a 2 , ...,a n ] 

where 

xi 1 — ^ Z, 

Xi 1 y BiZ — on for all 2 < i < n. 
Namely, for any polynomial h(x\,X2, ■ ■ • ,x n ) in fC[xi, X2, ■ ■ ■ , x n ], 

4>(h(xi,x 2 , ■ • ■ , x n )) = h(Z, B 2 Z - a>2, ■ ■ ■ , B n Z - a n ). 
The map is a ring isomorphism by virtue of the inverse map 

r\z) = x u 

4>^ 1 (ai) = BiX\ — Xi for all 2 < i < n, 

i.e., 4>~ 1 (h(Z, a 2 , ■ ■ ■ ,a n )) = h(xi,B 2 xi -x 2 ,-.-, B n xi - x n ). 

Since is a ring isomorphism, so is We can prove from gcd(/, g) = 1 that 

gcd(0(/),0( fl )) = l > i.e. 

gcd^/J.^lgcd^" 1 ^/))^- 1 ^))) = gcd(/, ff ) = 1. 

Now consider the Sylvester resultant 

p 1 (a 2 ,...,a n ) = Res z {(j)(f),(j)(g)) e K[a 2 , . . . , a n }. 

Because gcd(0(/), <p{g)) = 1, we have pi ^ 0. Remark that f(Z) = f(Z,B 2 Z—p2, ■ ■ ■ ,B n Z- 
p n ) and g{Z) = f(Z, B2Z — P2, ■ ■ . , B n Z — p n ) are obtained by substituting a, = pi for 
all 2 < i < n into </>(/) and (f>(g), respectively From the assumption that 

deg(f(Z,B 2 Z,...,B n Z)) = e 1 and deg(g(Z,B 2 Z, . . . ,B n Z)) = e 2 , 

it follows that dcg(/) = e\ and deg(g) = e 2 , and thus 

Res z (/,<?) =/Dl(p2,-..,Pn)- 

Therefore, gcd(/, g) = 1 is equivalent to pi(p2, ■ ■ ■ ,Pn) 0- The probability estimate then 
follows from Lemma 1 and the degree estimate dcg(pi) < 2 deg(/) deg(<?) = 2eie2- □ 

Now we are ready to present the main theorem. 

Theorem 2. Suppose that f,g 6 K,[xi_, . . . , x n ] with deg(/) = e\ < deg(g) = e 2 , and 
that / does not divide g, i.e., / j g. Suppose that B2, ■ ■ ■ , B n ,p 2 , . . . ,p n are chosen 
randomly and uniformly from a finite set W C /C. Let f,g £ be the univariate 

polynomials constructed from /, g as in (4). Then 

prob(/tp)> (i_^y (1 * y (i_ 

Proof. According to Lemma 2, the probability of picking B 2 , ■ ■ . , -B n randomly from 
such that deg(/) = ei and deg(g) = e2 is greater than (1 — yp^) • (1 — jf^r)- Suppose 
that f \ g, which is equivalent to that there exists one factor f\ of / with deg f\ > such 
that gcd(/i,<?) = 1. Remark that deg / = deg / implies that deg/i = deg/i. Having 
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gcd(/i,.g) = 1 and under the condition deg/i = deg/i and degg = degg, we easily get 
from Lemma 3 that 

2eie 2 



Prob(gcd(/i,a) ^ 1) = 1 - Prob(gcd(/ 1 ,.9) = 1) < 



\w\ 



If / has a factor /i such that gcd(/i, g) = 1 then clearly f \g. Therefore, the probability 
estimate of / { g is obtained as the product of two probabilities 

Prob({gcd(/,.g) = 1} | {deg/ = e x , degg = e 2 }) 

and Prob(deg/ = e\, degg = e 2 ). □ 

Remark 3. Stated in Theorem 2, our method is able to remove most of the non-invariant 
polynomials from the candidate invariants by applying univariate polynomial division, 
and the remaining non-invariant polynomials can be excluded by multivariate polynomial 
division. However, in practice, applying only univariate polynomial division can separate 
all the non-invariant polynomials from the actual invariants. 

In the end, let us analyze the complexity of our method of determining divisibility 
between multivariate polynomials based on linear transformation and univariate poly- 
nomial divisibility test. The extended Euclidean algorithm is usually applied to deter- 
mine divisibility between two univariate polynomials /, g G K,[x] with the complex- 
ity 0(deg(/) • deg(<7)). The complexity of combining linear transformation is then given 
by the following theorem. 

Theorem 3. Given /, g G /C[xi, . . . ,x n ] with deg(/) = e\ < deg(g) = e 2 . Assume that 
T = max(T(/), T(g)), where T(f) (resp. T{g)) denotes the number of terms in / (resp. g). 
Let / and g be the univariate polynomials as constructed in (4). Then the complexity of 
the method of determining divisibility between / and g, based on linear transformation 
and the extended Euclidean algorithm is 0(Xeg log 2 eg). In the worst case where both / 
and g are dense, the complexity is 0{{n + eg) 62 ). 

Proof. First, we analyze the cost for translating a multivariate polynomial / into the 
univariate polynomial f(Z) = f(Z, B 2 Z —p 2 , B„Z —p n ) as in (4). Let x^x^ 2 ■ ■ • xf™ 
be a term in the polynomial /, and denote e = X)"=i The expansion of Yi7=i(BiZ — 
Vi ) di can be c omputed by 0(e log 2 e) operations, according to the fan-in process method 



in iPanl (|2001l ). From the assumption, T is an upper bound of T(f) and T(g), and e 2 is 
an upper bound for the degrees of all of the terms in / and g. Therefore, the cost of 
computing f,g is bounded by 0(T e2 log 2 eg). If / and g are dense, we have e 2 <C T. 
In this case, the cost is bounded by 0((n + Z2Y 2 ) since T < (™^ 62 ) < (n + e 2 ) 62 , in 
which (™+ e2 ) is the maximum number of distinct terms of a polynomial in n variables 
and with a total degree bound e 2 . 

Clearly, the cost of applying the extended Euclidean algorithm on f(Z) and g(Z) is 
bounded by O(ejeg). In contrast to computing the linear transformation, this part of 
computation is negligible since e\ < T(f) < T in general. Therefore, the total cost of our 
method is 0(T eg log 2 eg), or is 0((n + eg) es ) if /, g are dense polynomials. □ 



Theorem 4. (jMonagan and Pearcd . 120071) The complexity for computing division of 



multivariate polynomials using heap is 0(NM log M), where M is the number of terms 
in the quotient and N is the number of terms in the divisor. 
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Remark 4. Suppose that r]\, . . . ,rj r are the candidates of polynomial invariants and there 
are t polynomials whose associated univariate polynomials satisfy rji(Z)\fji(Z'). In terms 
of the probability analysis shown in Theorem 2, all these t polynomials are invariants of 
the given program with high probability. In practice, the number of invariants is usually 
much less than that of the candidates obtained from the vanishing ideals of sample 
points, i.e. t <^ r. Therefore, compared with r times multivariate polynomial divisions, 
our method just needs r times univariate polynomial divisions, and t times multivariate 
polynomial divisions. Based on the complexity analysis in Theorems 3 and 4, our method 
of combining univariate polynomial division with partial multivariate polynomial division 
will be much more efficient. 

3.2. Generating Invariants of Polynomial Loops with Initial Values 

Based on the results in Section 3.1, we now present how to generate polynomial equa- 
tion invariants for polynomial loops with initial values. 

Let P be a polynomial loop program with n program variables, and e an upper bound 
for the total degree of its potential polynomial invariants. Then we need at most (™^ e ) 
sample points to determine a polynomial invariant in n variables with total degree 
bound e. We can obtain a finite set S containing no more than (™^ e ) sample points 
by executing the program P. 

Then, we apply Buchberger-Moller algorithm to compute the vanishing ideal 

= (Vi, ■•■,Vr) 

of the point set S, where rji, . . . ,r) r is a Grobner basis of the ideal I(S). In addition, we 
can get the minimal degree e'(< e) of all the polynomials in I(S). 

Consequently, we verify whether each candidate in {771 = 0,...,r] r = 0} is really 
an invariant. Clearly, all the r\i satisfy Initiation condition, since they belong to the 
vanishing ideal of sample points of the program. Therefore, the remaining task is to 
determine whether rji satisfies polynomial-scale consecution, i.e., r]i(V)\rn(V'). The 
techniques in Section 3.1 can be applied to carrying on this test efficiently. 

As a result, we can either generate the polynomial invariants of total degree < e or 
conclude that polynomial invariants with degree < e' do not exist. 

4. Algorithms 

We now present an algorithm to generate polynomial invariants of polynomial loops 
with initial values, and several examples are given to illustrate the algorithm. 

4-1. Polynomial Invariant Generation of Loop Programs with Numerical Initial Values 

By omitting the guard conditions of loop programs, the following algorithm states how 
to generate polynomial invariants of polynomial loops with numerical initial values. 
Algorithm InvGen 
Input: 

P : a polynomial loop program with numerical intial values 
n : the number of program variables 

e : an upper bound for the degree of polynomial invariants of P 
a : a term ordering on monomials of IC[xi, . . . , x n ] 
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Step 1: Set InvGen := NULL; 

Step 2: Translate the loop program P into a transition system T = (V, L, T, lo, 6); 
Step 3: Get a point set S containing (™^ e ) sample points by executing the loop P; 
Step 4: Apply Buchberger-Moller algorithm to compute a Grobner basis {r/i, 772, . . . , rj r } 
of the vanishing ideal I(S) of S with respect to the term ordering a. Then 

{ r/i = 0, . . . , ?7 r = } 

is a set of candidate invariants. Set e' = min{deg(77i), . . . , deg(77 r )}; 
Step 5: For i = 1, ...,r, verify whether 77^ = is really an invariant using Theo- 
rem 1: for r)i(V) and ^j(V') construct their corresponding univariate polynomials fji(Z) 
and fji(Z') as in (4); 

5.1 if fji(Z) \ fji(Z') then remove rji = from the candidates; 

5.2 otherwise if fji(Z)\fji(Z') then carry on multivariate polynomial division to check 
if Vi(V)\Vi(V). 

5.2.1 if r]i(V) \ r)i(V) then remove rji = from the candidates; 

5.2.2 otherwise if f?i(V)|77i(V) then InvGen := InvGen A {77, = 0}. 
Step 6: Return the polynomial equation invariants of the program P. 

6.1 If InvGen = NULL, then return "the polynomial equation invariants with degree < 
e' do not exist." 

6.2 Otherwise, return "the polynomial equation invariants of program P is InvGen." 
The following theorem is given to analyze the complexity of Algorithm InvGen. 

Theorem 5. The complexity of Algorithm- InvGen is 0((n + e)^ e+2 ), where n is the 
number of program variables and e is the upper bound for the degree of polynomial 
invariants. 

Proof. Here, we ignore the cost of Step 3, i.e., the time for executing the loop to obtain 
the point set S. In Step 4, we need to apply Buchberger-Moller algorithm to compute 
a vanishing ideal of (™^ e ) sample points, and by Remark 1, the complexity involved in 
Step 4 is 0(n 2 (n + e)4 e ). According to Theorem 3, the worst-case complexity in Step 5, 
where all the polynomials 771, . . . , rj r are dense, is bounded by 0(r e (n + e) 2e log(n + e)). 
Hence, the total cost of Algorithm- InvGen is bounded by 0((n + e)^ e+2 ). □ 

We give some working examples to illustrate Algorithm InvGen. 



Example 1 ((Petter, 2004 ; Rodn'guez-Carbonell and Kapur . 20071) ). Generate invari 



ants of the following program PI: 

(x,y) :=(0,0); 
lo : while ? do 

(x,y) := (x + y 5 ,y + 1) 
end while 

Set e = 7 to be the degree bound of the polynomial invariants of PI, and let a be the 
graded lexicographical ordering y -< x. 
Step 1 Set InvGen := NULL; 
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Step 2 We represent the program PI as a transition system (V, L, T, j^o}, ©) : 
{ 

V = {x,y} 
L = {lo} 

T = {t} with t = (lo, Iq, x' = x + y 5 A y' = y + 1, true) 
= {a; = OAy = O} 

} 

Step 3 By running the program PI with the initial values (0,0), we get (™^ e ) = 36 
sample points: 

(0,0), (0,1), (1,2), (33,3), (235306401,34), (280741825,35). 

Step 4 Apply Buchberger-Mollcr algorithm to get a vanishing ideal (rji, . . . ,r)6) of the 
above 36 sample points with respect to the graded lexicographical ordering y -< x. We 
only list 

r h (x,y) = -12x + 2y 6 -6y 5 + 5y i -y 2 , (5) 

since the other 5 polynomials are too huge. And, the minimal degree of polynomials 
in the vanishing ideal is e' = 6; 
Step 5 Determine if rji(x, y) \ r]i(x' , y') for i = 1, . . . , 6, and we get that 

rn{x, y) = -12x + 2y 6 - 6y 5 + 5y 4 - y 2 = 

is an invariant of the program PI, while the other 5 candidates are not invariants of 
the program PI. Moreover, we can conclude that the minimal degree of the polynomial 
invariants r](x, y) = is e' = 6. 

Now let us consider the loop programs where both the (numerical) initial values and 
guard conditions are taken into account. Then the candidate invariants can be verified by 
solving a semi-algebraic system. Indeed, if rji(V) = satisfies the consecution condition 
of inductive invariants, then 

Vi(V) = 0Ag T (V) 1=^0 = 0, 

i.e., each real solution of r]i(V) = A g T (V) also satisfies i]i(V') = 0, or equivalently, the 
semi-algebraic system 

(r h {V) = 0)Ag T (V)A( m (V')^0) 
has no real solutions, which can be verified by Maple package RegularChains. 

Accordingly, to generate polynomial invariants of a polynomial loop P with numerical 
initial values and guard condition, Step 5 of Algorithm- InveGen needs to be revised as 
follows: 

Step 5' For i = 1, . . . , r, verify whether the candidate rji = is an invariant by solving 
the semi-algebraic system 

(r H (V) = 0)Ag T (V)A(m(V')^0). 
If the above system has no real roots, then InvGen := InvGen A {?/; = 0}. 
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4-2. Polynomial invariant generation of loop programs with symbolic initial values 

Algorithm InvGen presents how to generate polynomial invariants of loop programs 
with numerical initial values. However, it is quite often in the program analysis and 
verification that only symbolic initial values of the program variables are given. Then 
the boolean value of the guard conditions will rely on the (concrete) initial values of 
program variables. In such situations, it is hard to compute exact sample points and 
thus Algorithm InvGen cannot be applied directly. 

Here, we supply two alternative ways to generate invariants for loop programs with 
symbolic initial values. One is to generate (symbolic) sample points by omitting the guard 
conditions and apply Algorithm InvGen with these sample points, as illustrated by the 
following example. 



Example 2 ( Kovacsl ( 2007 ). Example 5.35). Generate invariants of the program P2: 

(x,r) :=(f,0); 
while x > r do 

(x, r) := (x — r, r + 1); 
end while 

Unlike in Example 1, the initial value of x is a symbolic value, but the numerical 
comparison in the guard condition x > r relies on the concrete value of x. Instead, we 
omit the guard condition x > r. A set of sample points can be collected, for example, S = 
{(§ ,0), (§, 1), (§ -1,2), (f -3,3), (f -6,4)}. Applying Algorithm InvGen with a degree 
bound e = 2 of the polynomial invariants, we obtain an invariant 



2x + r 2 -r - a = 



of the program P2. 



An alternative approach is to combine Algorithm InvGen and rational function inter- 
polation. The idea is the following. 

Let Xi, . . . ,x n be the program variables of a polynomial loop P and iti, . . . , u m the 
symbolic initial values of P. Suppose that 

t] =p 1 (ui,...,u m )Ti H \-pk(ui,...,u m )T k = 

is an invariant of the program P, where Ti are monomials in xi,...,x n and pi are 
polynomials in u\, . . . ,u m - Actually, the representation of r\ is not unique, for example, 
for any ce/C\{0},c-?y = 0is also an invariant. To make the representation of r\ unique, 
we suppose that 

V = T 1 + ^T 2 + ■■■ + ^T k . 
Pi Pi 

The polynomial 77 can be obtained using Algorithm InvGen and rational function in- 
terpolation method as follows. First, we assign the initial values ui, . . . , u m to be random 
values. Then for each evaluation of the symbolic initial values u\,... ,u m , apply Algo- 
rithm InvGen to obtain a polynomial invariant which involves only xi,. .. At last, 
the polynomial r\ can be recovered by rational function interpolation based on variable 
by variable interpolation a nd early termination techniques. More details can be found in 
Kaltofen and Yand (|2007th 



Let us look at two examples to illustrate the above method. 
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Example 3 (jRodrfguez-Carbonell and Kapurl (|2007l ). Example 18). Generate invariants 
of the program P3: 

(x,y,u,v) := (a,b,b,a); 
while x ^ y do 
if x > y 

(x,y,u,v) 
else 

(x,y,u,v) 
end if 
end while 



(x -y,y,u,u + v); 
(x,y - x,u + v,v); 



In the program P3, we need to consider not only the guard condition x ^ y, but also 
two branch conditions x > y and x < y. Let (a, b) be assigned randomly, for example, 

287 751 
(ai ' 6l) = ( 253'890 } - 
An invariant for the initial values (ai, b\) can be generated using Algorithm InvGen: 

. 1 H2585 H2585 

mx, y,u,v) = l xu yv = 0. 

/v ' y ' ' ' 215537 215537 ,y 

Assume that the invariants to be found have the form: 

Vi K x , y,u,v,a,b) = l-\ — — xu + — - — —yv = 

Pi(a,b) Pi(a,b) 

where x,y,u,v are variables and a, b are parameters. 

Next, we recover the rational functions and . By evaluating (a, b) to the 

following numerical values, we get the corresponding invariants r)i for i = 2, 3, . . . : 



a-2 = 


93 
122 ' 


b 2 


_ 301 
992' 


m = 


i - 


903 x 903 yu u ' 


a.3 = 


349 
247' 


b 3 


239 
378' 


m = 


l - 


46683 46683,,,, _ n 
83411^" 83411 i> u u ' 


(24 = 


301 
6 ' 


in 


= 3, 


Vi = 


l - 


MT 2 ^ _ Wiy v = ^' 


a 5 = 


283 
352' 


b. 


17 

744 ' 


V5 


l - 


130944 130944 n 
4811 4811 » y — u ' 



By use of multivariate rational fu nction interpolation based on early termination tech- 
nique in (jKaltofen and YansJ . 120071) , an invariant of program P3 can be obtained: 

yv 



1 - 



2ab 2ab 



0, 



or, eq uivalently, — 2ab + xu + yv = 0. 



In (jRodriguez-Carbonell and Kapuii 120071 ). the authors obtained the same invariant 
as above using the fixed-point procedure, but they ignored the branch conditions x > y 
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and x < y. Our method takes into account the branch conditions, and are therefore more 
accurate in general. 



Example 4. By modifying the algorithm (jPetterl . 120041 ) for computing sums of powers, 
we consider the following more complicated series of programs with symbolic initial 
values: 

(x,y) := (a,b); 
while true do 

(x,y) := (x + y k ,y + 1); 
end while 

for k = 1,2,..., 15. We apply rational function interpolation method and Algorithm- 
InveGen to compute the polynomial invariants 

Vk(x,y,a,b) = 

that correspond to the (polynomial) assignments 

(x,y) := {x + y k ,y + l), 

for fc = l,2,...,15, respectively. The expressions of rjk are given in Table 1. 



5. Conclusions 

In this paper, we present a new method for generating polynomial equation invariants 
for polynomial loop programs. We first generate vanishing ideals of program sample points 
to get candidate invariants, then provide a probabilistic method to falsify divisibility 
so that we can exclude quickly non-invariants in a basis of the computed vanishing 
ideals. Our approach avoids first-order quantifier elimination and cylindrical algebraic 
decomposition as well as they do not depend on any abstraction interpretation methods. 
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V, 


a,b) 


= y 2 -y-2x + b-b 2 + 2a 


mix, 


li- 


a,b) 


= -6a; + y-3y 2 + 2y 3 -b + 3b 2 -2 b 3 + 6a 


mix, 


lt- 


a,b) 


= 4x + y 2 - 2y 3 + y 4 - b 2 + 2b 3 - b 4 +4a 


T)4,(X, 


li- 


a,b) 


= 6 y 5 - 30 x - y + 10 y 3 - 15 y 4 + b - 10 b 3 + 15 b 4 - 6 b 5 + 30 a 


Vb(x, 


lt- 


a,b) 


= -12 x- y 2 + 5y 4 - 6y 5 + 2y G + b 2 - 5b 4 + 6b 5 - 2b 6 + 12 a 


T]s(x, 


ti- 


a,b) 


= 21 y 5 - 42x + y + 6y 7 - 7 y 3 - 21y 6 - b+ 7b 3 - 21b 5 
+21 b 6 - 6 b 7 + 42 a 


rjr(x, 


ll- 


a, b) 


= -24x + 2y 2 -7y 4 + 14y 6 - 12y 7 + 3y 8 -2b 2 + 7b 4 
-14 b 6 + 12 b 7 - 3 6 s + 24 a 


mix, 


ti- 


a,b) 


= -90 x - 3 y + 20 y 3 - 42 y 5 + 10 y 9 + 60 y 7 - 45 y 8 + 3 b 
-20 b 3 + 42 6 5 - 60 b 7 + 45 b 8 - 10 b 9 + 90 a 




ll- 


a,b) 


= -3y 2 + 10y 4 - 14/ - 10y 9 + 2y 10 + 15y 8 + 3b 2 
-10b 4 + 14 b 6 - 15 6 8 + 10 b 9 - 2 6 10 + 20a 


Vw(x 


11 


a,b) 


= -66 x + 5 y - 33 y z + 66 y 5 + 55 y 9 - 33 y 10 - 66 y 7 + 6 y 11 
-5 6 + 33 b 3 - 66 b 5 + 66 b 7 - 55 b 9 + 33 6 10 - 6 b 11 + 66 a 


r/ii(x 


y 


a,b) 


= -24 x + 10 y 2 - 33 y 4 + 44 y 6 + 22 y w - 33 y 8 - 12 y 11 + 2 y 12 
-10b 2 +33 6 4 - 44 b 6 + 33b 8 - 22 b 10 + 12 6 11 - 2 b 12 +24 a 




u 


a,b) 


= -2730 1 - 691 y + 210 y 13 + 4550 y 3 - 9009 y 5 - 5005 y 9 
+8580 y + 2730 y 11 - 1365 y 12 + 691 6 - 4550 b 3 + 9009 b 5 
-8580 b 7 + 5005 b 9 - 2730 b 11 + 1365 b 12 - 210 b 13 + 2730 a 




u 


a,b) 


= -420 x - 210 y 13 - 691 y 2 + 2275 y 4 + 30 y 14 - 3003 y 6 
-1001 y 10 +2145 y 8 + 455 y 12 + 691 b 2 - 2275 b 4 + 3003 b 6 
-2145 b 8 + 1001 b 10 - 455 b 12 + 210 b 13 - 30 b 14 + 420 a 




u 


a,b) 


= -90 x + 105 y + 105 y 13 - 691 y 3 + 6 y 15 + 1365 y 5 + 715 y 9 

-45 y 14 - 1287 y 7 - 273 y 11 - 105 b + 691 b 3 - 1365 b 5 + 1287 b 7 
-715 b 9 + 273 b 11 - 105 b 13 + 45 b 14 - 6 b 15 + 90 a 


Vi5(x 


y 


a,b) 


= -48 x + 420 y 2 - 24 y 15 - 1382 y 4 + 60 y 14 + 1820 y 6 + 572 y 10 
-1287 y 8 +3y 16 - 182 y 12 - 420 b 2 + 1382 b 4 - 1820 b 6 + 1287 b 8 
-572 b 10 + 182 b 12 - 60 b 14 + 24 b 15 - 3 b ie + 48 a 



Table 1: The Invariants in Example 4 
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